#include <kuai/mol/atom.h>
#include <kuai/mol/bond.h>

#ifndef KUAIMISC_MOL_H
#define KUAIMISC_MOL_H

namespace kuai {

	class Molecule
		: Noncopyable
	{
	public:
		Molecule();

		explicit
			Molecule(
				const AtomArray& atoms
			);
		explicit
			Molecule(
				const AtomArray& atoms, const BondArray& bonds
			);

		void
			clear();
		void
			swap(
				Molecule& v0
			);

		const String
			formula() const;
		RealNumber
			weight() const;

		inline bool
			empty() const
		{
			return atoms.empty();
		}

		// The following functions were derived from Molecule Interface.
		Index
			countAtoms() const;
		Index
			countBonds() const;
		const Atom&
			getAtom(
				Index index
			) const;
		Atom&
			getAtom(
				Index index
			);
		const Atom&
			getAtom(
				const Bond& bond, Index index
			) const;
		Atom&
			getAtom(
				const Bond& bond, Index index
			);
		const Bond&
			getBond(
				Index index
			) const;
		Bond&
			getBond(
				Index index
			);
		const Bond*
			getBond(
				const Atom& a1, const Atom& a2
			) const;
		Bond*
			getBond(
				const Atom& a1, const Atom& a2
			);
		Index
			degree(
				const Atom& a
			) const;
		const Atom&
			getNeighborAtom(
				const Atom& ai, Index j
			) const;
		Atom&
			getNeighborAtom(
				const Atom& ai, Index j
			);
		const Bond&
			getNeighborBond(
				const Atom& a1, Index i
			) const;
		Bond&
			getNeighborBond(
				const Atom& a1, Index i
			);

		Index
			indexOf(
				const Atom& v0
			) const;
		Index
			indexOf(
				const Bond& v0
			) const;

	public:
		HashMap<String, String>								properties;

	private:
		AtomArray											atoms;
		BondArray											bonds;
		std::vector<std::vector<std::pair<Index, Index> > >	adjs;
		void setup();
	};

	typedef SharedPtr<Molecule>				MoleculePtr;

	template<typename MoleculeType>
	class SubMolecule
		: Noncopyable
	{
	public:
		typedef MoleculeType MolType;
		typedef typename MolType::AtomType AtomType;
		typedef typename MolType::BondType BondType;

	public:
		explicit
			SubMolecule(
				MolType& mol
			)
		{
			Index nAtoms = mol.countAtoms();
			atoms.reserve(nAtoms);
			for (Index i = 0; i < nAtoms; ++i) {
				atoms.push_back(&mol.getAtom(i));
			}
			Index nBonds = mol.countBonds();
			bonds.reserve(nBonds);
			for (Index i = 0; i < nBonds; ++i) {
				bonds.push_back(&mol.getBond(i));
			}
			setup(mol);
		};
		explicit
			SubMolecule(
				MolType& mol,
				const std::vector<AtomType*>& atoms0,
				const std::vector<BondType*>& bonds0,
				bool check_bond
			)
			: atoms(atoms0)
		{
			if (check_bond) {
				std::set<Atom*> filter(atoms0.begin(), atoms0.end());
				bonds.reserve(bonds0.size());
				for (std::vector<BondType*>::const_iterator
					it = bonds0.begin(); it != bonds0.end(); ++it)
				{
					Atom* atom1 = &mol.getAtom(**it, 0);
					Atom* atom2 = &mol.getAtom(**it, 1);
					if (filter.find(atom1) != filter.end()
						&& filter.find(atom2) != filter.end())
					{
						bonds.push_back(*it);
					}
				}
			}
			else {
				bonds = bonds0;
			}
			setup(mol);
		}

		explicit
			SubMolecule(
				MolType& mol,
				const std::vector<AtomType*>& atoms0
			)
			: atoms(atoms0)
		{
			Index nBonds = mol.countBonds();
			for (Index i = 0; i < nBonds; ++i) {
				BondType& bondI = mol.getBond(i);
				const AtomType& a1 = mol.getAtom(bondI, 0);
				const AtomType& a2 = mol.getAtom(bondI, 1);

				if (indexOf(a1) < atoms.size() && indexOf(a2) < atoms.size()) {
					bonds.push_back(&bondI);
				}
			}
			setup(mol);
		}
		explicit
			SubMolecule(
				MolType& mol,
				const std::vector<BondType*>& bonds0
			)
			: bonds(bonds0)
		{
			std::set<Atom*> atomsSet;
			for (std::vector<BondType*>::const_iterator
				i = bonds0.begin(); i != bonds0.end(); ++i)
			{
				AtomType& atom1 = mol.getAtom(**i, 0);
				AtomType& atom2 = mol.getAtom(**i, 1);
				atomsSet.insert(&atom1);
				atomsSet.insert(&atom2);
			}
			atoms.insert(atoms.end(), atomsSet.begin(), atomsSet.end());
			setup(mol);
		}

		explicit
			SubMolecule(
				SubMolecule& mol,
				const std::vector<AtomType*>& atoms0,
				const std::vector<BondType*>& bonds0,
				bool check_bond
			)
			: atoms(atoms0)
		{
			if (check_bond) {
				std::set<Atom*> filter(atoms0.begin(), atoms0.end());
				bonds.reserve(bonds0.size());
				for (std::vector<BondType*>::const_iterator
					it = bonds0.begin(); it != bonds0.end(); ++it)
				{
					Atom* atom1 = &mol.getAtom(**it, 0);
					Atom* atom2 = &mol.getAtom(**it, 1);
					if (filter.find(atom1) != filter.end()
						&& filter.find(atom2) != filter.end())
					{
						bonds.push_back(*it);
					}
				}
			}
			else {
				bonds = bonds0;
			};
			setup(mol);
		}

		explicit
			SubMolecule(
				SubMolecule& mol,
				const std::vector<AtomType*>& atoms0
			)
			: atoms(atoms0)
		{
			Index nBonds = mol.countBonds();
			for (Index i = 0; i < nBonds; ++i) {
				BondType& bondI = mol.getBond(i);
				const AtomType& a1 = mol.getAtom(bondI, 0);
				const AtomType& a2 = mol.getAtom(bondI, 1);

				if (indexOf(a1) < atoms.size() && indexOf(a2) < atoms.size()) {
					bonds.push_back(&bondI);
				}
			}
			setup(mol);
		}

		explicit
			SubMolecule(
				SubMolecule& mol,
				const std::vector<BondType*>& bonds0
			)
			: bonds(bonds0)
		{
			std::set<Atom*> atomsSet;
			for (std::vector<BondType*>::const_iterator
				i = bonds0.begin(); i != bonds0.end(); ++i)
			{
				AtomType& atom1 = mol.getAtom(**i, 0);
				AtomType& atom2 = mol.getAtom(**i, 1);
				atomsSet.insert(&atom1);
				atomsSet.insert(&atom2);
			}
			atoms.insert(atoms.end(), atomsSet.begin(), atomsSet.end());
			setup(mol);
		}


	public:
		Index
			countAtoms() const
		{
			return atoms.size();
		}
		Index
			countBonds() const
		{
			return bonds.size();
		}
		AtomType&
			getAtom(
				Index index
			)
		{
			return *atoms[index];
		}
		const AtomType&
			getAtom(
				Index index
			) const
		{
			return *atoms[index];
		}
		AtomType&
			getAtom(
				const BondType& bond, Index index
			)
		{
			assert (index == 0 || index == 1);
			assert (indexOf(bond) < countBonds());
			Index iB = indexOf(bond);
			assert (iB < countBonds());
			return *ends[iB*2+index];
		}
		const AtomType&
			getAtom(
				const BondType& bond,
				Index index
			) const
		{
			assert (index == 0 || index == 1);
			assert (indexOf(bond) < countBonds());
			Index iB = indexOf(bond);
			assert (iB < countBonds());
			return *ends[iB*2+index];
		}
		BondType&
			getBond(
				Index index
			)
		{
			return *bonds[index];
		}
		const BondType&
			getBond(
				Index index
			) const
		{
			return *bonds[index];
		}
		BondType*
			getBond(
				const AtomType& a1, const AtomType& a2
			)
		{
			assert (indexOf(a1) < countAtoms() && indexOf(a2) < countAtoms());
			Index i1 = indexOf(a1);
			for (Index j = 0; j < adj[i1].size(); ++j) {
				if (atoms[adj[i1][j].first] == &a2) {
					return bonds[adj[i1][j].second];
				}
			}
			return NULL;
		}
		const BondType*
			getBond(
				const AtomType& a1, const AtomType& a2
			) const
		{
			assert (indexOf(a1) < countAtoms() && indexOf(a2) < countAtoms());
			Index i1 = indexOf(a1);
			for (Index j = 0; j < adj[i1].size(); ++j) {
				if (atoms[adj[i1][j].first] == &a2) {
					return bonds[adj[i1][j].second];
				}
			}
			return NULL;
		}
		Index
			degree(
				const AtomType& a
			) const
		{
			assert (indexOf(a) < countAtoms());
			Index i1 = indexOf(a);
			return adj[i1].size();
		};
		AtomType&
			getNeighborAtom(
				const AtomType& a, Index j)
		{
			assert (indexOf(a) < countAtoms());
			assert (j < degree(a));
			Index i1 = indexOf(a);
			return *atoms[adj[i1][j].first];
		}
		const AtomType&
			getNeighborAtom(
				const AtomType& a, Index j
			) const
		{
			assert (indexOf(a) < countAtoms());
			assert (j < degree(a));
			Index i1 = indexOf(a);
			return *atoms[adj[i1][j].first];
		}
		BondType&
			getNeighborBond(
				const AtomType& a, Index j
			)
		{
			assert (indexOf(a) < countAtoms());
			assert (j < degree(a));
			Index i1 = indexOf(a);
			return *bonds[adj[i1][j].second];
		}
		const BondType&
			getNeighborBond(
				const AtomType& a, Index j
			) const
		{
			assert (indexOf(a) < countAtoms());
			assert (j < degree(a));
			Index i1 = indexOf(a);
			return *bonds[adj[i1][j].second];
		}

		Index
			indexOf(
				const AtomType& v0
			) const
		{
			std::vector<AtomType*>::const_iterator
				it = std::find(atoms.begin(), atoms.end(), &v0);
			if (it != atoms.end()) {
				return it - atoms.begin();
			}
			return Index(-1);
		}
		Index
			indexOf(
				const BondType& v0
			) const
		{
			std::vector<BondType*>::const_iterator  it = std::find(bonds.begin(), bonds.end(), &v0);
			if (it != bonds.end()) {
				return it - bonds.begin();
			}
			return Index(-1);
		}

		inline bool
			isEmpty() const
		{
			return atoms.empty();
		}

	private:
		std::vector<AtomType*> atoms;
		std::vector<BondType*> bonds;
		std::vector<AtomType*> ends;
		std::vector<std::vector<std::pair<Index, Index> > > adj;

	private:
		void
			setup(
				MolType& mol
			)
		{
			adj.resize(atoms.size());
			ends.reserve(bonds.size()*2);
			for (Index i = 0; i < bonds.size(); ++i)
			{
				AtomType& atom1 = mol.getAtom(*bonds[i], 0);
				AtomType& atom2 = mol.getAtom(*bonds[i], 1);
				ends.push_back(&atom1);
				ends.push_back(&atom2);
				Index i1 = indexOf(atom1), i2 = indexOf(atom2);
				assert (i1 != INVALID_INDEX && i2 != INVALID_INDEX);
				adj[i1].push_back(std::make_pair(i2, i));
				adj[i2].push_back(std::make_pair(i1, i));
			}
		}
		void
			setup(
				SubMolecule& mol
			)
		{
			adj.resize(atoms.size());
			ends.reserve(bonds.size()*2);
			for (Index i = 0; i < bonds.size(); ++i)
			{
				AtomType& atom1 = mol.getAtom(*bonds[i], 0);
				AtomType& atom2 = mol.getAtom(*bonds[i], 1);
				ends.push_back(&atom1);
				ends.push_back(&atom2);
				Index i1 = indexOf(atom1), i2 = indexOf(atom2);
				assert (i1 != INVALID_INDEX && i2 != INVALID_INDEX);
				adj[i1].push_back(std::make_pair(i2, i));
				adj[i2].push_back(std::make_pair(i1, i));
			}
		}
	};

	typedef SharedPtr<SubMolecule<Molecule> >	SubMoleculePtr;

	MoleculePtr clone(SubMoleculePtr v0);

	enum VisitedFlag
	{
		WHITE_FLAG,
		GRAY_FLAG,
		BLACK_FLAG
	};

	template<typename MoleculeType>
		class BasicMolecularVisitor
	{
	public:
		typedef MoleculeType MolType;
		typedef typename MoleculeType::AtomType AtomType;
		typedef typename MoleculeType::BondType BondType;

		void
			start(MoleculeType& mol, const AtomType&)
		{ }
		void finish(MoleculeType& mol, const AtomType&)
		{ }
		void find(MoleculeType& mol, const AtomType& vFind, const AtomType& vFrom)
		{ }
		void back(MoleculeType& mol, const AtomType& vFind, const AtomType& vFrom)
		{ }
	};

	template<typename VisitorType>
	void
		dft(
			typename VisitorType::MolType& mol,
			VisitorType& visitor,
			VisitedFlag flags[],
			const typename VisitorType::AtomType& v1
		)
	{
		visitor.start(mol, v1);
		flags[mol.indexOf(v1)] = GRAY_FLAG;
		Index nDegree = mol.degree(v1);
		for (Index i = 0; i < nDegree; ++i) {
			const typename VisitorType::AtomType& neighbor = mol.getNeighborAtom(v1, i);
			Index iNeighbor = mol.indexOf(neighbor);
			if (flags[iNeighbor] == WHITE_FLAG) {
				visitor.find(mol, neighbor, v1);
				dft(mol, visitor, flags, neighbor);
			}
			else {
				visitor.back(mol, neighbor, v1);
			}
		}
		flags[mol.indexOf(v1)] = BLACK_FLAG;
		visitor.finish(mol, v1);
	}
	template<typename VisitorType>
	void
		dft(
			typename VisitorType::MolType& mol,
			VisitorType& visitor,
			VisitedFlag flags[]
		)
	{
		Index nAtoms = mol.countAtoms();
		for (Index i = 0; i < nAtoms; ++i) {
			if (flags[i] == WHITE_FLAG) {
				dft(mol, visitor, mol.getAtom(i));
			}
		}
	}
	template<typename VisitorType>
	void dft(typename VisitorType::MolType& mol, VisitorType& visitor)  {
		std::vector<VisitedFlag> flags(mol.countAtoms(), WHITE_FLAG);
		dft(mol, visitor, &flags[0]);
	};
	template<typename VisitorType>
	void dft(typename VisitorType::MolType& mol, VisitorType& visitor, const typename VisitorType::AtomType& v1)  {
		std::vector<VisitedFlag> flags(mol.countAtoms(), WHITE_FLAG);
		dft(mol, visitor, &flags[0], v1);
	}

	template<typename VisitorType>
	void bft(typename VisitorType::MolType& mol, VisitorType& visitor, VisitedFlag flags[], const typename VisitorType::AtomType* initQueue[], size_t nQueue)
	{
		std::queue<Index> q;
		for (size_t i = 0; i < nQueue; ++i) {
			Index iAtom = mol.indexOf(*initQueue[i]);
			q.push(iAtom);
			flags[iAtom] = GRAY_FLAG;
		}
		for (; !q.empty(); q.pop()) {
			Index i = q.front();
			const typename VisitorType::AtomType& a = mol.getAtom(i);
			visitor.start(mol, a);
			Index nDegree = mol.degree(a);
			for (Index j = 0; j < nDegree; ++j) {
				const typename VisitorType::AtomType& neighbor = mol.getNeighborAtom(a, j);
				Index iNeighbor = mol.indexOf(neighbor);
				if (flags[iNeighbor] == WHITE_FLAG) {
					visitor.find(mol, neighbor, a);
					q.push(iNeighbor);
					flags[iNeighbor] = GRAY_FLAG;
				}
				else {
					visitor.back(mol, neighbor, a);
				}
			}
			flags[i] = BLACK_FLAG;
			visitor.finish(mol, a);
		}
	}

	template<typename VisitorType>
	void bft(typename VisitorType::MolType& mol, VisitorType& visitor, VisitedFlag flags[], const typename VisitorType::AtomType& v1)  {
		const typename VisitorType::AtomType* atoms[] = { &v1 };
		bft(mol, visitor, flags, atoms, 1);
	}
	template<typename VisitorType>
	void bft(typename VisitorType::MolType& mol, VisitorType& visitor, const typename VisitorType::AtomType& v1)  {
		std::vector<VisitedFlag> flags(mol.countAtoms(), WHITE_FLAG);
		bft(mol, visitor, &flags[0], v1);
	}
	template<typename VisitorType>
	void bft(typename VisitorType::MolType& mol, VisitorType& visitor, VisitedFlag flags[])  {
		Index nAtoms = mol.countAtoms();
		for (Index i = 0; i < nAtoms; ++i) {
			if (flags[i] == WHITE_FLAG) {
				bft(mol, visitor, flags, mol.getAtom(i));
			}
		}
	}
	template<typename VisitorType>
	void bft(typename VisitorType::MolType& mol, VisitorType& visitor)  {
		std::vector<VisitedFlag> flags(mol.countAtoms(), WHITE_FLAG);
		bft(mol, visitor, &flags[0]);
	}

	template<typename MoleculeType>
	class ConnectionPartVisitor
		: public BasicMolecularVisitor<MoleculeType>
	{
	public:
		typedef SubMolecule<MoleculeType> ResultType;
		typedef SharedPtr<SubMolecule<MoleculeType> > ResultPointer;
		typedef std::vector<ResultPointer> ResultPointerArray;

	public:
		ConnectionPartVisitor(Index reserved) {
			atoms.reserve(reserved);
		}

		void find(MoleculeType& mol, const AtomType& vFind, const AtomType& vFrom)
		{
			atoms.push_back(const_cast<AtomType*>(&vFind));
		}

		void clear() {
			atoms.clear();
		}

		ResultPointer getResult(MoleculeType& mol) {
			return ResultPointer(new SubMolecule<MoleculeType>(mol, atoms));
		};

		const std::vector<AtomType*>& getResult() const {
			return atoms;
		};

	public:
		static Index countParts(MoleculeType& mol) {
			Index nAtoms = mol.countAtoms();
			std::vector<VisitedFlag> flags(nAtoms, WHITE_FLAG);
			Index result = 0;
			ConnectionPartVisitor<MoleculeType> visitor(nAtoms);
			for (Index i = 0; i < nAtoms; ++i) {
				if (flags[i] == WHITE_FLAG) {
					visitor.clear();
					dft(mol, visitor, &flags[0], mol.getAtom(i));
					++result;
				}
			}
			return result;
		}

		static ResultPointerArray getConnectedParts(MoleculeType& mol) {
			Index nAtoms = mol.countAtoms();
			std::vector<VisitedFlag> flags(nAtoms, WHITE_FLAG);
			ResultPointerArray result;
			ConnectionPartVisitor<MoleculeType> visitor(nAtoms);
			for (Index i = 0; i < nAtoms; ++i) {
				if (flags[i] == WHITE_FLAG) {
					visitor.clear();
					visitor.atoms.push_back(&mol.getAtom(i));
					dft(mol, visitor, &flags[0], mol.getAtom(i));
					result.push_back(visitor.getResult(mol));
				}
			}
			return result;
		}

		static ResultPointerArray getConnectedParts(SubMolecule<MolType>& mol) {
			Index nAtoms = mol.countAtoms();
			std::vector<VisitedFlag> flags(nAtoms, WHITE_FLAG);
			ResultPointerArray result;
			ConnectionPartVisitor<SubMolecule<MoleculeType> > visitor(nAtoms);
			for (Index i = 0; i < nAtoms; ++i) {
				if (flags[i] == WHITE_FLAG) {
					visitor.clear();
					visitor.atoms.push_back(&mol.getAtom(i));
					dft(mol, visitor, &flags[0], mol.getAtom(i));
					result.push_back(ResultPointer(new ResultType(mol, visitor.getResult())));
				}
			}
			return result;
		}

	public:
		std::vector<AtomType*> atoms;
	};

	template<typename MoleculeType>
	class SpanTreeVisitor
		: public BasicMolecularVisitor<MoleculeType>
	{
	public:
		virtual void find(MoleculeType& mol, const AtomType& vFind, const AtomType& vFrom) {
			sources[&vFind] = &vFrom;
		}

	public:
		const AtomType& getSource(const AtomType& v) const {
			assert (sources.find(&v) != sources.end());
			return *(sources.find(&v)->second);
		}

	protected:
		HashMap<const AtomType*, const AtomType*> sources;
	};

	template<typename MoleculeType>
	class RingFinder
		: public SpanTreeVisitor<MoleculeType>
	{
	public:
		typedef SubMolecule<MoleculeType> ResultType;
		typedef SharedPtr<ResultType> ResultPointer;
		typedef std::vector<ResultPointer> ResultPointerArray;

	public:
		enum {
			UNDECIDED = INVALID_INDEX,	// Maybe in a huge ring (large than max rank) or not on a ring.
			NOT_RING = 0				// Not on a ring.
		};

	public:
		static Index findMinRing(MoleculeType& mol, const BondType& b, std::vector<AtomType*>& result, Index maxRank0)
		{
			AtomType& atom1 = mol.getAtom(b, 0);
			AtomType& atom2 = mol.getAtom(b, 1);
			RingFinder visitor(atom2, atom1, maxRank0);
			std::vector<VisitedFlag> flags(mol.countAtoms(), WHITE_FLAG);
			flags[mol.indexOf(atom1)] = flags[mol.indexOf(atom2)] = GRAY_FLAG;
			try {
				bft(mol, visitor, &flags[0], atom2);
				return NOT_RING;
			}
			catch (RingFinder& finder) {
				if (finder.finded()) {
					finder.getResult(result);
					return result.size();
				}
				else {
					return Index(UNDECIDED);
				}
			}
		}

		static ResultPointerArray getRingCandidates(MoleculeType& mol)
		{
			SubMolecule<MoleculeType> sub(mol);
			return getRingCandidates(sub);
		}

		static ResultPointerArray getRingCandidates(SubMolecule<MoleculeType>& mol)
		{
			Index nAtoms0;
			ResultPointer candidates = ResultPointer(new ResultType(mol));
			do {
				nAtoms0 = candidates->countAtoms();
				std::vector<MoleculeType::AtomType*>  atoms;
				atoms.reserve(nAtoms0);
				for (Index i = 0; i < nAtoms0; ++i) {
					MoleculeType::AtomType* atomI = &(candidates->getAtom(i));
					Index d = candidates->degree(*atomI);
					if (d > 1) {
						atoms.push_back(atomI);
					}
				}
				candidates = ResultPointer(new ResultType(*candidates, atoms));
			} while (candidates->countAtoms() < nAtoms0);

			ResultPointerArray results = ConnectionPartVisitor<MoleculeType>::getConnectedParts(*candidates);
			PredNotRingCandidate pred;
			ResultPointerArray::iterator it = std::remove_if(results.begin(), results.end(), pred);
			results.erase(it, results.end());
			return results;
		}

		static void getRingRankOfBonds(MoleculeType& mol, Index result[], Index maxRank0) {
			ResultPointerArray candidates = getRingCandidates(mol);
			// memset(result, INVALID_INDEX, sizeof(Index)*mol.countBonds());
			memset(result, 0, sizeof(Index)*mol.countBonds());
			for (ResultPointerArray::const_iterator it = candidates.begin(); it != candidates.end(); ++it) {
				Index nAtoms = (*it)->countAtoms();
				Index nBonds = (*it)->countBonds();
				if (nAtoms == nBonds) {
					for (Index j = 0; j < nBonds; ++j) {
						BondType& bondJ = (*it)->getBond(j);
						assert (mol.indexOf(bondJ) != INVALID_INDEX);
						result[mol.indexOf(bondJ)] = nBonds;
					}
				}
				else {
					assert (nAtoms < nBonds);
					std::vector<Index> resultI;
					getRingRankOfBondsOfCandidates(**it, resultI, maxRank0);
					for (Index j = 0; j < nBonds; ++j) {
						Index bondJ = mol.indexOf((*it)->getBond(j));
						result[bondJ] = resultI[j];
					}
				}
			}
		}

		static ResultPointerArray getRingBlocks(MoleculeType& mol, Index maxRank0) {
			ResultPointerArray candidates = getRingCandidates(mol);
			ResultPointerArray results; results.reserve(candidates);

			for (ResultPointerArray::const_iterator it = candidates.begin(); it != candidates.end(); ++it) {
				Index nAtoms = it->countAtoms();
				Index nBonds = it->countBonds();
				if (nAtoms == nBonds) {
					if (nAtoms <= maxRank0) {
						results.push_back(*it);
					}
				}
				else {
					assert (nAtoms < nBonds);
					std::vector<Index> resultI;
					getRingRankOfBondsOfCandidates(*it, resultI, maxRank0);
					std::vector<BondType*> bonds;
					for (size_t j = 0; j < resultI.size(); ++j)
					{
						if (3 <= resultI[j] && resultI[j] <= maxRank0) {
							bonds.push_back(&it->getBond(j));
						}
					}
					sub = new SubMolecule<MolType>(*it, bonds);

					ResultPointerArray blocks = ConnectPartVisitor<MolType>::getConnectedParts(sub);
					results.insert(results.end(), blocks.begin(), blocks.end());
				}
			}
			return results;
		}

		static ResultPointerArray getRingBlocks(SubMolecule<MoleculeType>& mol, Index maxRank0) {
			ResultPointerArray candidates = getRingCandidates(mol);
			ResultPointerArray results; results.reserve(candidates.size());

			for (ResultPointerArray::iterator it = candidates.begin(); it != candidates.end(); ++it) {
				Index nAtoms = (*it)->countAtoms();
				Index nBonds = (*it)->countBonds();
				if (nAtoms == nBonds) {
					if (nAtoms <= maxRank0) {
						results.push_back(*it);
					}
				}
				else {
					assert (nAtoms < nBonds);
					std::vector<Index> resultI;
					getRingRankOfBondsOfCandidates(**it, resultI, maxRank0);
					std::vector<BondType*> bonds;
					for (size_t j = 0; j < resultI.size(); ++j)
					{
						if (3 <= resultI[j] && resultI[j] <= maxRank0) {
							bonds.push_back(&((*it)->getBond(j)));
						}
					}
					SubMolecule<MolType> sub(**it, bonds);

					ResultPointerArray blocks = ConnectionPartVisitor<MolType>::getConnectedParts(sub);
					results.insert(results.end(), blocks.begin(), blocks.end());
				}
			}
			return results;
		}

		static ResultPointerArray getSingleRings(MoleculeType& mol, Index maxRank0) {
			ResultPointerArray blocks = getRingBlocks(mol, maxRank0);
			ResultPointerArray result; result.reserve(blocks.size());
			for (ResultPointerArray::const_iterator it = blocks.begin(); it != blocks.end(); ++it) {
				if ((*it)->countAtoms() == (*it)->countBonds()) {
					result.push_back(*it);
				}
				else {
					assert ((*it)->countAtoms() < (*it)->countBonds());
					ResultPointerArray rings = getSingleRingsFromRingBlock(**it, maxRank0);
					result.insert(result.end(), rings.begin(), rings.end());
				}
			}

			return results;
		}


		static ResultPointerArray getSingleRingsFromRingBlock(SubMolecule<MoleculeType>& mol, Index maxRank0) {
			ResultPointerArray result;
			Index nBonds = mol.countBonds();

			std::vector<Index> ranks;
			ranks.resize(mol.countBonds(), Index(UNDECIDED));
			std::vector<ResultType::AtomType*> atoms;
			std::vector<ResultType::BondType*> bonds;
			for (Index i = 0; i < nBonds; ++i) {
				Bond& bondI = mol.getBond(i);
				if (ranks[i] == UNDECIDED) {
					Index nRank = RingFinder<ResultType>::findMinRing(mol, bondI, atoms, maxRank0);
					if (3 <= nRank && nRank <= maxRank0) {
						bonds.clear();
						bonds.reserve(nRank);
						bonds.push_back(&bondI);
						if (nRank < ranks[i]) {
							ranks[i] = nRank;
						}
						for (Index j = 1; j < nRank; ++j) {
							Bond* bondJ = mol.getBond(*atoms[j-1], *atoms[j]);
							bonds.push_back(bondJ);
							assert (bondJ != NULL);
							Index idJ = mol.indexOf(*bondJ);
							if (ranks[idJ] > nRank) {
								ranks[idJ] = nRank;
							}
						}
					}
					assert (atoms.size() == bonds.size());
					result.push_back(ResultPointer(new SubMolecule<MoleculeType>(mol, atoms, bonds, false)));
				}
			}
			return result;
		}

	private:
		class PredNotRingCandidate {
		public:
			inline bool operator()(const ResultPointer& frag) const
			{
				Index nAtoms = frag->countAtoms();
				Index nBonds = frag->countBonds();
				return nBonds < nAtoms;
			}
		};

		static inline bool isChainBond(ResultType& mol, typename ResultType::BondType& bond) {
			typename ResultType::AtomType& atom1 = mol.getAtom(bond, 0);
			typename ResultType::AtomType& atom2 = mol.getAtom(bond, 1);
			return mol.degree(atom1) <= 2 || mol.degree(atom2) <= 2;
		}

		static void getRingRankOfBondsOfCandidates(ResultType& mol, std::vector<Index>& result, Index maxRank0) {
			Index nBonds = mol.countBonds();
			result.resize(nBonds, Index(UNDECIDED));
			std::vector<ResultType::AtomType*> atoms;
			for (Index i = 0; i < nBonds; ++i) {
				BondType& bondI = mol.getBond(i);
				if (result[i] == UNDECIDED || !isChainBond(mol, bondI)) {
					Index nRank = RingFinder<ResultType>::findMinRing(mol, bondI, atoms, maxRank0);
					if (3 <= nRank && nRank <= maxRank0) {
						if (nRank < result[i]) {
							result[i] = nRank;
						}
						for (Index j = 1; j < nRank; ++j) {
							BondType* bondJ = mol.getBond(*atoms[j-1], *atoms[j]);
							assert (bondJ != NULL);
							Index idJ = mol.indexOf(*bondJ);
							if (result[idJ] > nRank) {
								result[idJ] = nRank;
							}
						}
					}
				}
			}
		}

	private:
		RingFinder(const AtomType& a1, const AtomType& a2, Index maxRank0) {
			assert (maxRank0  >= 3);
			source = const_cast<AtomType*>(&a1);
			target = const_cast<AtomType*>(&a2);
			result = NULL;
			maxRank = maxRank0-2;
			depth.insert(std::make_pair(source, 0));
			depth.insert(std::make_pair(target, 0));
		}

	public:

		void back(MoleculeType& mol, const AtomType& vFound, const AtomType& vFrom)
		{
			SpanTreeVisitor::back(mol, vFound, vFrom);
			if (target == &vFound && &vFrom != source) {
				result = const_cast<AtomType*>(&vFrom);
				throw *this;
			}
		}

		virtual void find(MoleculeType& mol, const AtomType& vFind, const AtomType& vFrom)
		{
			SpanTreeVisitor::find(mol, vFind, vFrom);
			assert (depth.find(const_cast<AtomType*>(&vFind)) == depth.end());
			assert (depth.find(const_cast<AtomType*>(&vFrom)) != depth.end());
			Index depthI = depth.find(const_cast<AtomType*>(&vFrom))->second + 1;
			depth.insert(std::make_pair(const_cast<AtomType*>(&vFind), depthI));
			if (depthI > maxRank) {
				throw *this;
			}
		}

		void getResult(std::vector<AtomType*>& atoms) const {
			assert (result != NULL);
			atoms.clear();
			atoms.reserve(10);
			atoms.push_back(target);
			for (AtomType* i = result; i != source; i = const_cast<AtomType*>(&getSource(*i))) {
				atoms.push_back(i);
			}
			atoms.push_back(source);
		}

		bool finded() const {
			return result != NULL;
		}

	private:
		AtomType* source;
		AtomType* target;
		AtomType* result;
		Index			maxRank;
		HashMap<AtomType*, Index> depth;
	};
}

namespace std {
	template<>
	inline void
		swap<kuai::Molecule>(
			kuai::Molecule& v1, kuai::Molecule& v2
		)
	{
		v1.swap(v2);
	}

}

#endif
